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ABSTRACT 

We show how to develop an expansion of nearly oblate systems in terms of a set of 
potential-density pairs. A harmonic (multipole) structure is imposed on the potential 
set at infinity, and the density can be made everywhere regular. We concentrate on 
a set whose zeroth order functions describe the perfect oblate spheroid of de Zeeuw 
(1985). This set is not bi-orthogonal, but it can be shown to be complete in a weak 
sense. Poisson's equation can be solved approximately by truncating the expansion 
of the potential in such a set. A simple example of a potential which is not one of 
the basis functions is expanded using the symmetric members of the basis set up to 
fourth order. The basis functions up to first order are reconstructed approximately 
using 10,000 particles to show that this set could be used as part of an TV-body code. 

Key words: methods: analytical, numerical - celestial mechanics, stellar dynamics - 
galaxies: kinematics and dynamics 



INTRODUCTION 
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o 

in i 

\jr^ The use of harmonic basis function expansion of spher- 
^■^ical or nearly spherical systems has been described by 
Clutton-Brock (1973), Polyachenko & Shukhman (1981) 
<™H and, amongst others, Hernquist & Ostriker (1992). Such an 
^""expansion is a useful tool for finding approximate solutions 



Q to Poisson's equation 
P, 



(1) 



where p is the density and $ the potential. They can be used 
in simulations of collisionless, and/or collisional particles, in- 
cluding gas (Hernquist & Ostriker 1992, Barnes & Hernquist 
1993, Weil & Hernquist 1993) and in normal mode analy- 
sis (Robijn 1994). Solutions have the advantage of being 
smooth and efficient to generate, especially if the expansion 
can be truncated at low order. 

The functions consist of a set of potential-density pairs 
{$jim, Pjim}, where (j, I, m) are integers. The potential and 
density of a general system may be written 



$(r) 



p(r) 



jlm 



y ] AjlmPjlr, 

jlm 



(r) 



(2) 



(3) 



yhere 



V $j im = Pjlrn , (4) 

By convention j labels the radial dependence of $ji m , and I 
and m label the polar dependence. In spherical coordinates 
typically 

$ jim = J F ji (r) J P im (cos0)e"^, (5) 



where Pi m is the associated Legendre polynomial. In 
the works mentioned above Fji can be chosen to make 
jlm, p jlm} bi-orthogonal, i.e 



where 



A 



jlm 



is the Kronecker delta, and hence 



$ jim (r) / 9(r)rf 3 



Let us write 

F j1 (t) = W j1 (t)U j1 (t), 



(6) 



(7) 



(8) 



where V ' 3 i — > constant and Wji has the desired asymptotic 
form of $ji m as r — > 0, oo,. We refer to a set in which 
Wji is chosen to have a multipole structure at infinity, as a 
'harmonic set,' viz. 



W jt (r) 



1 



Clutton-Brock (1973) chose 
zeroth-order potential and 

w n(r) = (1 + r2)i+1/2 ; 

Hernquist and Ostriker (1992) 



(1 + r 



a+i 



(9) 

a Plummer model for the 
(10) 

chose a Hernquist model and 
(11) 



The importance of the zeroth-order model lies in the ability 
to truncate the expansion at low order for systems in which 
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$(r) ~ $000. 

Saha (1993) points out that considerable flexibility may 
be gained by dropping the requirement of bi-orthogonality. 
In this case 



QikmPjlnd r = Sij klmn i 



(12) 



and the matrix Sijkimn must be inverted to find {Aji m }. 
Inversion of the matrix is eqivalent to Gram-Schmidt or- 
thogonalisation of the chosen function set. If the expansion 
is being truncated at low order, the matrix inversion is not 
costly, and only has to be done once in any case. The added 
flexibility of this approach arises in the freedom to spec- 
ify Woo arbitrarily, and thus lower order truncation can be 
done in systems which are still almost spherical, but are not 
nearly the same as a Plummer or a Hernquist model. 

In this paper we describe a class of function sets that 
have oblate or disk-like zeroth-order density. We concen- 
trate on a set which has a de Zeeuw (1985) perfect oblate 
spheroid as the zeroth-order model. We take the view that 
bi-orthogonality is a luxury, and find the inverse of Sijkimn 
in low order for the perfect oblate spheroid set. Finally we 
show that solution of Poisson's equation for a set of 10,000 
particles is feasible using the same set. 



2 A HARMONIC OBLATE SET 

We will use a prolate spheroidal polar coordinate system 
throughout, for which we define the coordinates (£, rj) in the 
(R, z) plane as follows: 



R 2 



e 2 (f 2 



(13) 



We see that £ £ [l,oo) is radius-like, and that r] £ [— 1, 1] is 
like cos 9 in spherical polar coordinates. Surfaces of constant 
£ are prolate spheroids with foci ati? = 0,2 = ±l, and they 
are asymptotically spherical of radius £ as £ — > oo. The 
foci (£ = 1,jj = ±1) are singular points of the coordinate 
system, which we may think of as the end-points of the 
'needle' £ = 1. 

In these coordinates the Laplacian can be written as 



e 2 V 2 =L(t, V ) + L( V ,0 + 



1 



d 2 



(t 2 - 1)(1 - V 2 ) d<j> 2 
where <f> is the usual polar variable in the x-y plane, and 



(14) 



HZ,*) 



1 



d 



(1-t) 



d 



(15) 



We recognise in £(£,»?) a close relative of the 9 part of V 2 
in spherical coordinates. A consequence of this similarity is 
that 



Qjlm(r) = Pjm(£)Plm(l])e 



irrwfi 



(16) 



is an eigenfunction of V 2 , where Pi m is the associated Leg- 
endre function. Unfortunately the corresponding density is 
not regular at the points of the needle. The eigenfunctions 
of V 2 which are regular are not expressible in closed form 
(Abramowicz and Stegun 1964). Thus we are motivated to 
seek a potential-density set which is regular, but not orthog- 
onal. 



One well-known potential in this coordinate system is 
de Zeeuw's perfect oblate spheroid: 



$pos 



£tan 1 (e£) — »;tan 1 (ejj) 



f 2 - V 2 

with the remarkably simple and regular density 



Ppos 



2e(l + e 2 



(1 + e 2 f 2 ) 2 (l + eV) 2 
or, in cylindrical coordinates 



Ppos 



Po 



(l + R 2 /(l + e 2 ) + z 2 



(17) 



(18) 



(19) 



(de Zeeuw 1985, Binney and Tremaine 1987). 

As a generalization of (17) let us consider potentials of 
the form 



h(t,y)-h( V ,t) 

£ 2 — Vj 2 



(20) 



This, and the associated density, is regular at the points of 
the needle provided 

Ht,v) = h(-t,-v) (21) 

(see appendix). 

Thus we seek a set of potential functions based on 



h jlm (t,v) = Wj^^XjmiOYimiv), (22) 

where Xj m and Yi m combine to give us the 'completeness' 
that we are looking for in a potential basis set. We demand 
that Xoo = loo = 1, so that the function Woo determines 
the zeroth-order potential $o- This should be chosen to 
match the system we are modelling closely, in which case 
the function set expansion can be truncated at low order; 
here we concentrate on 



Woo (£,»?) =£tan- 1 (e£), 



(23) 



giving a perfect oblate spheroid. Following Hernquist and 
Ostriker (1992), we tailor the function Wji to give our set a 
harmonic structure at infinity. We require that 



1 

£i + l : 



as £ — > oo. 



(24) 



It will also be useful to require that Wji and at least its first 
two derivatives are regular everywhere (see appendix). 

The choice of the functions Wji, X jrn and Yi m will be 
influenced by the requirement that p = V 2 <£> is regular. The 
inclusion of the index m in equation (22) indicates that we 
intend to extend the (R, z) plane into three dimensions by 
multiplying $ji m by e' m ^ . Having done this we find that the 
<f> part of V 2 gives the simplest contribution to pji m , viz. 



1 



d 2 



(£ 2 - 1)(1 - v 2 ) dfr 



(e-m-v 2 ) 



jlm- 



(25) 



Thus we have to be careful about the 2-axis (£ = 1 and 
r\ = ±1), because V 2 has a vanishing denominator there 
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(as in the spherical case). The way we deal with this in the 
spherical case, where cos 9 has an i-like operator in V 2 , is to 
make $ proportional to P; m (cos 9). In the spheroidal case, 
both £ and r\ have the operator L in V 2 , so we set 



being careful to define Pi m (x) in such a way as to make it 
real for x > 1. (Most of the {Pi m (x)} are not analytic at 
x = 1, so analytic continuation makes no sense, but all we 
need is a function that ~ Pi m (x) as x — > 1.) An obvious 
choice is 



2 \m\ rf' m ' 

Pim(x) = a(m) |l - x 2 \ m , Pi(x), 

where Pi(x) is the usual Legendre polynomial and 



a(m) 



(-l) m , \x\<l 

1, otherwise. 



(26) 



(27) 



a(m) is chosen to give the derivatives of Pi m the same sign 
on either side of i = ±1. The appendix shows that this 
choice of X jrn and Yi m gives a regular density on the 2-axis. 

Having chosen X jrn and Yi m , we turn to Wji and 
to the points of the needle. We know that Pi m (x) = 
( — l) l+rn Pim( — x), and ho(S,,ri) has been chosen to be sym- 
metric, so to ensure the regularity of p 3 im at the points of 
the needle we require that 

w 3l (Ln) = (-iy +l w 3l (-t-n). (28) 

Assuming that Woo (£,,v) i s symmetric, an example would 
be 



(3 + 1 



(1 + eHe + V 2 )) 3 



-Woo{tv) 



(29) 



Putting all the above together, a well-behaved set whose 
zeroth order potential is the perfect oblate spheroid is given 
by 



hjim(Z, v) 



p + i (30) 
f —r- ftan-^ef) P jm (S,)Pim(v)- 

Other sets with regular density can be trivially derived from 
this by replacing £tan -1 (e£) with some other symmetric 
function of £ and/or rj. 

Unfortunately, none but the zeroth-order density func- 
tion (18) can be expressed in compact form. 



3 RESULTS 

The four lowest order density functions of the perfect oblate 
spheroid set are shown in Figures 1 and 2. The zeroth- order 
density is smooth and everywhere positive. The higher order 
functions are more oscillatory as one would expect, and are 
sometimes negative. All the density functions have finite 
or zero mass since they are regular and have a multipole 
structure at infinity. 

Table 1 lists the matrix elements of Sijkimn in a square 
grid for an eccentricity e of 0.5. Because of the orthogonality 



of the perfect oblate spheroid function set in the <f> direction, 
we only need consider m = n; and since $ji m = for m > 
min[j, I] we can ignore larger values of m. The map between 
the new index p and (jlm) is defined in the table. The 
resulting matrix is still quite sparse because Sijkimn = 
if k + I is odd, by symmetry. The integrations were all 
carried out over a finite volume bounded by the curve £ = 3e. 
All elements are accurate to approximately 10 significant 
figures. 

Consider the potential-density pair ($, p) in which 



$(r) = $ooo(ei,r) + a$ oo(e 2 ,r), 



(31) 



where ei,e2 are eccentricities, and a is a constant. We can 
try to approximate $(r), for instance, by expanding in the 
function set {$ji m (ei , r)}. Thus we construct the vector 



bik 



<&ifcm(ei,r)p(r), 



and find {A t k m } by solving the linear equations 



klmn In — Oik 



(32) 



(33) 



jln 



The result of doing this in the case (ei,e2) = (0.5,1.0), 
a = 0.1 is shown in Figure 3. The values of Aji m were found 
for the symmetric (I even, m = 0) functions up to order 
(j) = (4,4) (15 functions in all). The solid curve shows 
the density p along the 2-axis. The lower dashed curve is 
the zeroth order approximation; the upper dashed curve is 
the expansion using the values (jlm) = {000, 020, 040, 100}; 
and the expansion using all 15 functions is indistinguishable 
from the solid curve. 

Consider a swarm of particles with positions {r a } drawn 
from a distribution proportional to \p q \, and with masses 
{m a }, all equal in magnitude, but negative wherever the 
density is negative. The quantity 



bp = ^ m a $ p = S p 



(34) 



should be approximately equal to the matrix element S pq 
and the quantity 



q-1 O 



(35) 



should be approximately equal to the Kronecker delta. In 
Table 2 we show the results of testing this assertion for the 
lowest five density functions using 10,000 particles. It shows 
that 8p q is approximately equal to 8 pq , at least within the 
expected tolerance, given that statistical errors are expected 
to be a few percent. We might think of this as a minimum 
requirement of feasibility for the function set to be used in 
particle simulations. 



4 DISCUSSION 

An aesthetic drawback of using function sets which are not 
bi-orthogonal is that we are sometimes not sure quite how 
complete they are, and we hesitate to call our set a 'ba- 
sis set' until we know what space it spans. Of course one 
can always make the tautological statement that {$jim} is 
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Figure l.The four lowest order densities of the set defined by equation (30) with e = 0.5 in the i j plane: a) ( . . ••• j = (0, 0, 0) is the 
perfect oblate ellipsoid; b) (0,1,0); c) (1,0,0); and d) (1,1,0). 



a complete basis of the subspace of all functions spanned 
by {$jim}- The question is then, is this subspace a use- 
ful one; can it represent a sufficiently broad class of sensible 
potential-density pairs for practical purposes? The spherical 
function sets referred to in Section 1 are the eigenfunctions 
of a hermitian operator, hence they are bi-orthogonal, and 
they form a complete bi-normal set (i.e. any function for 
which $V 2 $ is integrable can be represented by a linear 
combination of the set). For practical purposes formal com- 



pleteness is not an issue of great urgency, for one can only 
ever use a finite subset of even a rigorously complete basis. 
The best one can do might be, as in the present work, to 
present a set of functions that closely resembles some zeroth 
order model of interest, in the hope of finding the smallest 
possible subspace which is practically useful. We present 
some arguments below which justify the choice of hji m in 
equation (30) by appealing to a rather weak formal com- 
pleteness. 
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p 



p 




Figure 2. As Figure 1 with e = 10 in the (x,z) plane: a) (j, I, m) 
d) (1,1,0). 



Saha (1993) considers sets of the form (8), but in which. 
Wji does not depend on j (as is the case in all the examples 
mentioned). Following him, we consider the completeness 
of such sets from the point of view of approximation in the 
mean (Courant and Hilbert 1937, Chapter II). Consider the 
radial part of the Ith harmonic, and denote the error in the 
radial part of the potential by SFi(r). If Uji(r) is a poly- 
nomial in some variable u(r) with finite range and domain, 



P 



P 




(0,0,0) is the perfect oblate ellipsoid; b) (0,1,0); c) (1,0,0); and 



then it can be shown that the average error 

may be made arbitrarily small by increasing the order of 
the polynomial (i.e. by truncating the series at larger j). 
This proves completeness of {U } i} over the space of square 
integrable, piecewise continuous functions of u (ibid.). Thus 
{Uji} is complete in a very definite sense, but not neces- 
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Table 1. The matrix S made square and indexed according to p. Each element is accurate to approximately 10 significant figures 



(jlm) p 1 2 3 4 5 

(000) 1 0.45815368137 0.89588371437 

(100) 2 0.5186873810 1.283914361 

(010) 3 1.162423755 3.626567178 

(110) 4 1.313280882 3.754218875 

(111) 5 3.677108173 



0.3 

0.2 
P 
0.1 




Figure 3. The density along the z-axis of a simple model (solid 
line), and successive approximations to it using the perfect oblate 
spheroid set (dashed lines). 



Table 2. The matrix 8 calculated using 10,000 particles. 

(jlm) p 1 

(000) 1 1.005 

(100) 2 -0.019 

(010) 3 -0.050 

(110) 4 -0.018 

(111) 5 -0.057 



2 

-0.011 
1.029 
0.042 

-0.045 



3 

-0.001 
-0.0005 
1.01 



4 

0.004 
-0.011 



5 

0.0003 
-0.003 



-0.010 -0.005 



-0.003 1.019 0.006 



-0.013 0.006 



-0.005 



1.01 



sarily {fji}. For a finite mass distribution, the {fji} are 
necessarily harmonic and so completeness of {Uji} should 
be sufficient for practical purposes. 

Notice the denominator in (29) above. It was chosen 
so that in the spherical limit (e — > 0), $ji m tends to a set 
which is complete in the sense just described. The spherical 
limit is also approached at large radius by sets with finite e. 
Specifically, in this limit, 



jlm 



(1 + r 2 



tan 1 (r) 



1 + r 2 



Pi m (cos#)e"^.(37) 



This is in the form of (5) with Fji in the form of (8), and 
hence the weak completeness is ensured. 

In fact, we can go a stage further: writing Wji in the 
form of Clutton-Brock (1973) we would extract 



(38) 



We denote the corresponding sets of Clutton-Brock by 
{Fob} and {Ucb}- The {Ucb} turn out to be a particular 
set of ultraspherical, or Gegenbauer, polynomials. Clearly 
the {Uji} are of finite range and their domain is the same 
as the {Ucb}- The two sets are both linearly independent 
and for every member of {Uji} there is a member of {Ucb}- 
This is no more than a necessary condition for completeness 



of {Uji}, but completeness of {Uji} would imply complete- 
ness of {fji} because {Fcb} is complete. This gives us some 
confidence that {fji} might be of practical use. 

A simpler looking, but less successful choice of $ji m 
would be as the above, with (f 2 + »? 2 ) replacing the quantity 
in brackets in the denominator of h. This is similarly well- 
behaved in terms of its density, but in the limit e — > 



jlm 



(1 + r 2 



tan 1 (r) Pim(cos 8) e'^ , 



(39) 



in which we see that {$jim} are not even linearly indepen- 
dent (they do not depend on j), and the corresponding set 
of potential-density pairs thus falls far short of being com- 
plete. (It could only be used to model systems with a very 
restricted radial dependence.) 

In the spherical case the requirement that the poten- 
tial is harmonic at leads to a set in which each member is 
bi-normal and has finite mass. For the works mentioned in 
the Section 1 it also leads to a bi-orthogonal set. For the 
present purposes we would admit that the harmonic sruc- 
ture is not strictly necessary, but it has features which may 
be considered attractive. For example, each member of the 
potential set is again automatically bi-normal and has fi- 
nite mass. Also, members of the density set with I / all 
have zero mass as a consequence of their formal similarity 
to spherical multipoles at infinity. Thus, all the mass in a 
given distribution resides in the I = part of the expansion. 
Also, the higher order parts of the expansion bear a familiar 
relationship to their I = counterparts. Thus, for instance, 
our expectation that an I = 1 mode should be 'dumb bell' 
shaped, or dipolar, is confirmed in Figure 1. 

It is evident from Figure 2 that the members of the 
perfect oblate spheroid set with I ^ have a density which 
is considerably less flattened than those with 1 = 0. This 
problem is probably generic to the coordinate system, which 
becomes progressively less similar in shape to the zeroth 
order model as e increases. The perfect oblate ellipsoid set 
is unlikely to be practically useful in the limit of very flat 
systems. The density in all the members falls off as a power- 
law in radius. Thus they are also unlikely to be successful 
at reproducing realistic disc models, which have exponential 
decay of the density in the ^-direction. 

Many elliptical galaxies are thought to have density 
cusps in their inner parts. The spherical bi-orthogonal set of 
Hernquist and Ostriker (1992) is particularly successful at 
reproducing cuspy density distributions owing to the pres- 
ence of a singularity in the density set at the origin, associ- 
ated with the coordinate singularity there. All of the perfect 
oblate spheroid density set are smooth and continuous, and 
so are unlikely to be successful at reproducing cuspy density 
distributions. 

A mass distribution that has constant density p{r 2 ) on 
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the self similar ellipsoids 



constant 



2 2 

x y 
~ + T 2 



a 



has a potential given by 
du 



where 



T 2 (u) 



and 
A(«) 



A(«) 



p(t z )dt z 



>(u) 



y 



+ u b 2 + u 



+ u 



^/(a 2 + u)(b 2 + u)(c 2 + u) 



(40) 



(41) 



(42) 



(43) 



(Chandresekhar 1969, see also de Zeeuw 1985). So an at- 
tempt at constructing a basis set in ellipsoidal coordinates 
with a cusp at the origin could be made, along similar lines 
to the Hernquist model, by setting 

1 d f(r) 



P{r) 
Then 



t dr 1 + t 2 



A(«) 



(A + u)(n + u)(v + u) 



f(r(u)) du. 



(44) 



(45) 



(The perfect ellipsoid appears as the case f(r) = constant.) 
However, it is difficult to choose f(r) in such a way as to 
make (45) integrable, while retaining the singularity in p. 

Finally, we note that the a triaxial set of potential den- 
sity pairs could be constructed by a very similar method. 
Consider the ellipsoidal coordinate system (A, where 
the coordinates are the roots for u of the equation t 2 (u) = 
const (see for example de Zeeuw 1985). The coordinate sys- 
tem defined by equation (13) is the case 

2 7.2 

— = — = 1 + e 2 

c 2 c 2 

The perfect triaxial ellipsoid in these coordinates has a po- 
tential of the form 



pte 



F(X) 



(A 



-P)(A- 



F{y) 



(46) 



(n-v)(n-X) (y ■ 



so, by analogy with equation (20 
sider a potential set of the form 



\)( v -nY 

we would be lead to con- 



$ : 



h(X, fi, v) 
' (A - p)(A - v) 
h(p, A, v) 



h(v, A, fi) 



(47) 



(p — v)(ii — A) (v-\)(v-ii)' 



This form guarantees regularity of the potential at the sin- 
gular points of the coordinate system. 
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5 APPENDIX 

In this appendix we will show first the condition on h(^,r]) 
which makes $ regular at the points of the needle ({,»?) = 
(1,-1). Fet us expand V 2 $ about the points ({,»?) = 
(1, ±1) in a Taylor series in rj. In fact, regularity is trivially 
satisfied at r] = +1, so we only give the result for r] = — 1: 



Thus $ is regular at ({,»?) = (1,-1) provided 
h(l,-l) = h(-l,l). 

Next let us expand p = V 2 $ about (£, rj) = (1, -1): 



H 



2(1 -y) 2 
3# + 2(ftii(l,-l)-ftii(-l,l)) 



(Al) 



(A2) 



(A3) 



4(1 - y) 



+ 0(l-y) 



phere 



H = G + /j 10 (-l,l) + /iio(l,-l) 
- Aoi(-l,l) + /ioi(l,-l), 
and 

G = h(l,-l) - h(-l,l). 



(A4) 



(A5) 



In addition to (A2) then, the extra conditions for p to be 
regular are 

h 10 (-l, 1) + h 10 (l, -1) - /ioi(-l, 1) + h 01 (l, -1) = 0, (A6) 



/in(l,-l) - Aii(-l,l) 



(At) 



To satisfy all the regularity conditions (A2), (A6) and (A7) 
it is sufficient that h(^,r]) = h(—^,—r]). 



Next we show that X r 



Y 

+ 7 r 



make p regular on the 2-axis, £ = 1, r] < 1 and £ > 1, r] = ±1. 
First let us define 



-if 



(A8) 



And then we write 
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Kt,v) = W(t, V )P j m(t)P,m( V ), 



(A9) 



where at least the first two derivatives of FF(£, rj) are regular, 
and notice that 



-J^Qjrn(i) + (regular stuff) . 



(A10) 
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Now 



1 , r r , „ m 2 



V 2 $ = (Lt - L„) $ $, 

f 2 - »7 2 " j (« 2 -l)(l-»? 2 ) ' 

which we can write as 



1 1 
+ 



1 1 - i] 2 



(All) 



(e-m-v 2 ). 

+ (regular stuff) . 
The quantity in braces vanishes, and hence p is reguiar. 
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